##Analisis exploratorio

Consumo Dataset

## 'data.frame':    257 obs. of  27 variables:
##  $ X                 : int  3 4 5 6 7 8 9 10 11 12 ...
##  $ Anio              : int  2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 ...
##  $ Mes               : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ GasolinaSuper     : num  308157 307766 331910 315648 319668 ...
##  $ GasolinaRegular   : num  202645 205531 229500 210680 208164 ...
##  $ TotalGasolinas    : num  510802 513297 561410 526328 527832 ...
##  $ Diesel            : num  634667 642381 699807 586804 656948 ...
##  $ DieselLS          : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ DieselULS         : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ TotalDiesel       : num  634667 642381 699807 586804 656948 ...
##  $ GLP               : num  194410 174711 189234 174331 191745 ...
##  $ GasolinaAviacion  : num  1426 1458 1503 1561 1642 ...
##  $ Kerosina          : num  64026 62660 61362 61814 54098 ...
##  $ TurboJet          : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Bunker            : num  296767 328116 368590 396300 449369 ...
##  $ Asfalto           : num  48446 50597 27593 53794 60137 ...
##  $ PetCoke           : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ AceitesLubricantes: num  NA NA NA NA NA NA NA NA NA NA ...
##  $ GrasasLubricantes : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Solventes         : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Naftas            : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Ceras             : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ CrudoNacional     : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Butano            : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Orimulsion        : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ MezclasOleosas    : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Total             : num  1750545 1773220 1909499 1800933 1941772 ...

Podemos ver que todas las variables tratadas en las datasets son variables cuantitativas, por lo que haremos histogramas de las variables.

## consumo$Anio : 
##         Frequency Percent Cum. percent
## 2020           12     4.7          4.7
## 2019           12     4.7          9.3
## 2018           12     4.7         14.0
## 2017           12     4.7         18.7
## 2016           12     4.7         23.3
## 2015           12     4.7         28.0
## 2014           12     4.7         32.7
## 2013           12     4.7         37.4
## 2012           12     4.7         42.0
## 2011           12     4.7         46.7
## 2010           12     4.7         51.4
## 2009           12     4.7         56.0
## 2008           12     4.7         60.7
## 2007           12     4.7         65.4
## 2006           12     4.7         70.0
## 2005           12     4.7         74.7
## 2004           12     4.7         79.4
## 2003           12     4.7         84.0
## 2002           12     4.7         88.7
## 2001           12     4.7         93.4
## 2000           12     4.7         98.1
## 2021            5     1.9        100.0
##   Total       257   100.0        100.0

Variable X: no sigue distribución normal, y es lo esperado ya que X es una variable diferente para los datos, como un ID.

Variable Anio: No sigue una distribucion normal tampoco, pero nos muestra que hay una frecuencia mayor en los primeros años comenzando en 2000, luego las frecuencias en los demas años son similares hasta 2021, cuando los datos son muy pocos.

Variable mes: No sigue una distribución normal tampoco.

Variable Gasolina Super: No sigue una distribución normal.

Variable Gasolina Regular: No sigue una distribución normal.

Variable Gasolina Diesel: Si sigue una distribución normal, por lo que pueda que tenga un comportamiento peculiar en el estudio de las series lineales.

####Cruce de variables Al estar trabajando con las variables de gasolina, en los cruces de variables los haremos siempre con las variables de gasolinaSuper, Regular y Diesel.

De este cruce de variables podemos ver que puede que exista un tipo de tendencia, que con los años, el consumo de gasolina super creció con los años.

Podemos ver que lo mismo sucedió con gasolina Regular y los años, pero otra vez, es la diesel la que tiene un comportamiento irregular.

Ahora vamos a hacer el cruce de variables pero con meses.

plot(y=consumo$GasolinaSuper,x=consumo$Mes,ylab="Gasolina Super", xlab="Mes")

plot(y=consumo$GasolinaRegular,x=consumo$Mes,ylab="Gasolina Regular", xlab="Mes")

plot(y=consumo$TotalDiesel,x=consumo$Mes,ylab="Gasolina Diesel", xlab="Mes")

Ahora vemos que las variables se comportan de una manera aleatoria conforme a los meses del año. Ya sea cualquier tipo de las gasolinas tratadas.

Analisis Importacion

Breve resumen del data set de importaciones

## [1] "Dimensiones"
## [1] 245  29
## [1] "Resumen"
##        X              Anio           Mes         GasolinaSuper    
##  Min.   :  3.0   Min.   :2001   Min.   : 1.000   Min.   : 170293  
##  1st Qu.: 64.0   1st Qu.:2006   1st Qu.: 3.000   1st Qu.: 357350  
##  Median :261.0   Median :2011   Median : 6.000   Median : 435416  
##  Mean   :267.5   Mean   :2011   Mean   : 6.429   Mean   : 464691  
##  3rd Qu.:433.0   3rd Qu.:2016   3rd Qu.: 9.000   3rd Qu.: 558403  
##  Max.   :623.0   Max.   :2021   Max.   :12.000   Max.   :1227174  
##                                                                   
##  GasolinaRegular  TotalGasolinas        Diesel           DieselLS      
##  Min.   : 81015   Min.   : 274500   Min.   : 229765   Min.   : 691066  
##  1st Qu.:194830   1st Qu.: 558276   1st Qu.: 620918   1st Qu.: 891341  
##  Median :282509   Median : 718114   Median : 767285   Median :1056569  
##  Mean   :341982   Mean   : 806673   Mean   : 782290   Mean   :1056457  
##  3rd Qu.:469437   3rd Qu.:1021878   3rd Qu.: 904976   3rd Qu.:1195728  
##  Max.   :896841   Max.   :1861582   Max.   :1595699   Max.   :1592580  
##                                     NA's   :41        NA's   :204      
##    DieselULS      TotalDiesel           GLP         GasolinaAviacion 
##  Min.   : 1203   Min.   : 229765   Min.   :100562   Min.   :   48.0  
##  1st Qu.: 7176   1st Qu.: 655715   1st Qu.:210871   1st Qu.:  487.2  
##  Median :17701   Median : 800223   Median :371098   Median : 2132.0  
##  Mean   :21905   Mean   : 829959   Mean   :375817   Mean   : 3292.6  
##  3rd Qu.:38300   3rd Qu.: 994338   3rd Qu.:510507   3rd Qu.: 3367.5  
##  Max.   :48946   Max.   :1630636   Max.   :960841   Max.   :27979.1  
##  NA's   :225                                        NA's   :89       
##     Kerosina         TurboJet          Bunker           Asfalto        
##  Min.   :  1000   Min.   : 18373   Min.   :   8485   Min.   :   74.88  
##  1st Qu.: 33762   1st Qu.: 45061   1st Qu.: 172130   1st Qu.: 2996.07  
##  Median : 51860   Median : 69728   Median : 294609   Median : 5497.36  
##  Mean   : 54898   Mean   : 73330   Mean   : 308321   Mean   : 7231.09  
##  3rd Qu.: 67657   3rd Qu.: 95020   3rd Qu.: 417863   3rd Qu.:10237.67  
##  Max.   :184094   Max.   :158719   Max.   :1051764   Max.   :48364.50  
##  NA's   :50       NA's   :183                        NA's   :10        
##     PetCoke       AceitesLubricantes GrasasLubricantes   Solventes    
##  Min.   :  6993   Min.   :15086      Min.   :115.8     Min.   : 3207  
##  1st Qu.:149221   1st Qu.:18020      1st Qu.:313.8     1st Qu.: 6264  
##  Median :165003   Median :21592      Median :388.2     Median : 8881  
##  Mean   :218887   Mean   :22826      Mean   :436.4     Mean   : 9974  
##  3rd Qu.:244990   3rd Qu.:24292      3rd Qu.:546.0     3rd Qu.:11978  
##  Max.   :849463   Max.   :44746      Max.   :996.5     Max.   :24501  
##  NA's   :127      NA's   :216        NA's   :216       NA's   :216    
##      Naftas            Ceras           Butano       PetroleoReconstit
##  Min.   :  0.380   Min.   :107.9   Min.   :  0.37   Min.   :356364   
##  1st Qu.:  3.263   1st Qu.:297.5   1st Qu.: 72.29   1st Qu.:360619   
##  Median : 54.000   Median :404.4   Median : 72.29   Median :367573   
##  Mean   :138.722   Mean   :428.8   Mean   : 76.75   Mean   :489360   
##  3rd Qu.:280.017   3rd Qu.:522.2   3rd Qu.: 81.30   3rd Qu.:718114   
##  Max.   :350.920   Max.   :920.1   Max.   :149.09   Max.   :730957   
##  NA's   :233       NA's   :218     NA's   :206      NA's   :225      
##       MTBE         Orimulsion     MezclasOleosas   PetroleoCrudo
##  Min.   : 8184   Min.   :249983   Min.   : 166.7   Min.   :191  
##  1st Qu.:11062   1st Qu.:312030   1st Qu.: 966.5   1st Qu.:191  
##  Median :12680   Median :315916   Median :1668.6   Median :191  
##  Mean   :12952   Mean   :315469   Mean   :1626.9   Mean   :191  
##  3rd Qu.:14090   3rd Qu.:325083   3rd Qu.:2417.1   3rd Qu.:191  
##  Max.   :19431   Max.   :344685   Max.   :3134.8   Max.   :191  
##  NA's   :234     NA's   :232      NA's   :223      NA's   :244  
##   TotalMensual    
##  Min.   :1381787  
##  1st Qu.:2065597  
##  Median :2468706  
##  Mean   :2558891  
##  3rd Qu.:2924855  
##  Max.   :4322010  
## 
## [1] "Estructura"
## 'data.frame':    245 obs. of  29 variables:
##  $ X                 : int  3 4 5 6 7 8 9 10 11 12 ...
##  $ Anio              : int  2001 2001 2001 2001 2001 2001 2001 2001 2001 2001 ...
##  $ Mes               : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ GasolinaSuper     : num  373964 243091 312084 285055 300914 ...
##  $ GasolinaRegular   : num  177777 123116 161726 127339 168730 ...
##  $ TotalGasolinas    : num  551740 366207 473811 412394 469644 ...
##  $ Diesel            : num  566102 489526 575560 437745 552609 ...
##  $ DieselLS          : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ DieselULS         : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ TotalDiesel       : num  566102 489526 575560 437745 552609 ...
##  $ GLP               : num  194066 170703 161837 163049 171519 ...
##  $ GasolinaAviacion  : num  820 3054 677 3399 585 ...
##  $ Kerosina          : num  33834 67440 31787 25801 45529 ...
##  $ TurboJet          : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Bunker            : num  214582 294609 315264 205653 278371 ...
##  $ Asfalto           : num  27749 7504 26304 7886 8443 ...
##  $ PetCoke           : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ AceitesLubricantes: num  NA NA NA NA NA NA NA NA NA NA ...
##  $ GrasasLubricantes : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Solventes         : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Naftas            : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Ceras             : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Butano            : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ PetroleoReconstit : int  715344 370166 360530 359527 723346 360369 360570 719303 360635 717717 ...
##  $ MTBE              : int  8402 NA NA 8184 12680 NA 10642 13357 NA 19431 ...
##  $ Orimulsion        : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ MezclasOleosas    : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ PetroleoCrudo     : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ TotalMensual      : num  2312639 1769209 1945770 1623638 2262727 ...

Aqui se encuentran los histogramas de las variables cuantitativas para ver si hay normalidad aun que sea en las continuas. Podemos observar que las variables: GasolinaSuper, TotalGasolinas, Diesel y TotalDiesel siguen una distribución normal.

Cruze de variables con el anio y mes En las tres gráficas anteriores del cruce de variables de las importaciones de cualquier tipo de gasolina y el año, se puede observar que hay una tendencia a crecer. Esto implica que conforme pasa el tiempo, cada vez se importa más gasolina sin importar su tipo.

Por otro lado, al observar las gráficas del cruce de variables de cualquier tipo de gasolina con los meses, no se puede observar alguna tendencia. Esto implica que sin importar el mes, siempre va a haber una importación de combustible más o menos constante.

Para comenzar con el analisis de series de tiempo nos encontramos con la necesidad de unir la columna de mes y anio en una sola columna en formato de fecha, luego crear diferentes data frames para los 3 tipos de gasolinas para hacer una linea del tiempo de cada una.

####Linea de tiempo para consumo de gasolina super Primero se convierte el data frame en una serie de tiempo. Luego podemos observar el inicio, fin, frecuencia y vista de la linea del tiempo.

## [1] 2000    1
## [1] 2021    5
## [1] 12

De esta linea podemos ver a primera vista que posee una tendencia a crecer, y que hay momentos en la linea de tiempo que hace que las varianzas tambien se vean afectadas y sean diferentes entre si. Podemos ver que se trata de una serie con frecuencia anual donde hay datos de todos los meses desde enero de 2000 hasta mayo del 2021.Como se puede ver la cantidad de galones vendidos va aumentando a medida que pasan los años.

#####Descomposicion de la serie

Podemos observar una serie con tendencia a aumentar, que no es estacionaria en varianza, y además tiene estacionalidad.

Ahora vemos los datos de inicio y fin del train y del test

## [1] 2000    1
## [1] 2014   12
## [1] 2015    1
## [1] 2021    5

#####Estimar los parámetros del modelo. Cómo no es estacionaria en varianza le haremos una transformación logaritmica para hacerla constante en varianza.

Veremos si logramos estacionarizar la serie con esta transformación:

Podemos observar que las varianzas si se pudieron hacer mas constantes, sin embargo, se pueden observar picos de varianza desiguales, lo cual puede ser un problema y concreta que no se pudo estacionarizar en varianza en su totalidad.

Por lo que ahora trataremos de quitar la tendencia y estacionarizar en normalidad.

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: 0.3673
##   P VALUE:
##     0.7319 
## 
## Description:
##  Fri Aug 06 21:42:55 2021 by user: jfdel
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     DF: 0.3673
##   P VALUE:
##     t: 0.7895 
##     n: 0.773 
## 
## Description:
##  Fri Aug 06 21:42:55 2021 by user: jfdel

Ambas pruebas de Dickey-Fuller Test nos demuestra que p es mayor a 0.05 por lo que no podemos rechazar Ho, y por ende, asumir que no hay raices unitarias. Por lo que la serie de tiempo no es estacionaria en media.

Por lo que debemos de hacer lo mismo pero con una diferenciacion en la serie de tiempo.

## Warning in adfTest(diff(train)): p-value smaller than printed p-value
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: -17.2389
##   P VALUE:
##     0.01 
## 
## Description:
##  Fri Aug 06 21:42:55 2021 by user: jfdel

Como se puede ver el valor de p ahora sí está por debajo de 0.05 por lo que se puede descartar la hipótesis nula de que existan raices unitarias.

Solo se hizo uso de una diferenciación, por lo que d es 1, pero para conocer p y q se debe de realizar un gráfico de autocorrelación y autocorrelación parcial.

acf(logGS,50) 

pacf(logGS,50) 

Del ACF podemos ver que la gráfica no se anula y tiene un decrecimiento rapido. pero en el PACF se anula luego de 1, por lo que q es 1. P no se anula en ningun momento. por lo que proponemos un modelo, podemos pensar que p es 4.

Por lo que podemos proponer un modelo ARIMA con p =4 q=1 d = 1

Veamos si hay estacionalidad en la serie. Si volvemos a ver la descomposición de la serie en su componente estacional.

Al parecer sí existe estacionalidad en la serie. Para tener una idea de los parámetros estacionales veremos las funciones de autocorrelación y autocorrelación parcial con 36 resagos para ver en que momentos son los picos estacionales. Se usará la serie estacionarizada.

Cómo podemos ver en la función de autocorrelación tenemos un decaimiento en los picos estacionales, 12, 18, 24, 30, 36, y en la función de autocorrelación parcial, un pico significativo en el retardo 12. Eso sugiere que los parámetros del componente estacional son P = 2, D = 1 (el período es de 6 meses por lo que se puede hacer una diferenciación estacional), y Q=1.

Creamos los analisis del modelo Hacemos un analisis de coeficientes del AutoArima y Arima.

## Warning: package 'lmtest' was built under R version 4.0.5
## 
## Attaching package: 'lmtest'
## The following object is masked from 'package:epiDisplay':
## 
##     lrtest
## 
## z test of coefficients:
## 
##       Estimate Std. Error  z value  Pr(>|z|)    
## ar1  -0.530999   0.255698  -2.0767 0.0378324 *  
## ma1  -0.016694   0.240298  -0.0695 0.9446149    
## ma2  -0.462657   0.139363  -3.3198 0.0009008 ***
## sar1 -0.021098   0.080054  -0.2635 0.7921275    
## sar2 -0.137446   0.080886  -1.6993 0.0892688 .  
## sma1 -0.999970   0.089828 -11.1321 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## z test of coefficients:
## 
##       Estimate Std. Error  z value  Pr(>|z|)    
## ma1  -0.694455   0.061480 -11.2957 < 2.2e-16 ***
## sma1  0.403711   0.079137   5.1014 3.371e-07 ***
## sma2  0.218524   0.067869   3.2198  0.001283 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Vemos que todos los coeficientes son significativos y vamos a hacer un analisis de residuos.

Análisis de Residuales Los residuales deben ser lo más parecidos posible a un ruido blanco. Esto es que sean aleatorios (no están correlacionados entre sí), y que estén distribuidos normalmente. Analizamos los residuos de nuestro modelo fitArima

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,2)(2,1,1)[12]
## Q* = 13.882, df = 18, p-value = 0.7367
## 
## Model df: 6.   Total lags used: 24

Según los gráficos podemos ver que la distribución de los datos parece ser normal, y que no hay correlaciones significativas. Según el test de Ljung-Box los datos se distribuyen de forma independiente puesto que el p-value es mayor a 0.05, por lo que no se puede rechazar la hipótesis nula. Esto significa que el modelo generado es aceptable para predecir.

Analizando los residuos del modelo generado de forma automática por R tenemos los siguientes resultados:

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,1)(0,0,2)[12]
## Q* = 35.203, df = 21, p-value = 0.02682
## 
## Model df: 3.   Total lags used: 24

El modelo autogenerado no es un modelo aceptable para predecir.

##Predicción con el modelo generado

##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Jan 2015       13.03072 12.96888 13.09256 12.93614 13.12530
## Feb 2015       13.00463 12.93677 13.07249 12.90085 13.10841
## Mar 2015       13.09947 13.02944 13.16950 12.99237 13.20658
## Apr 2015       13.04365 12.96995 13.11735 12.93094 13.15636
## May 2015       13.02189 12.94554 13.09824 12.90512 13.13866
## Jun 2015       12.97006 12.89072 13.04940 12.84872 13.09139
## Jul 2015       13.04983 12.96784 13.13183 12.92443 13.17523
## Aug 2015       13.05805 12.97337 13.14273 12.92854 13.18755
## Sep 2015       13.01030 12.92308 13.09752 12.87690 13.14370
## Oct 2015       13.04912 12.95939 13.13885 12.91189 13.18635
## Nov 2015       13.02783 12.93569 13.11997 12.88691 13.16875
## Dec 2015       13.16829 13.07378 13.26281 13.02374 13.31284
## Jan 2016       13.04320 12.94569 13.14071 12.89407 13.19233
## Feb 2016       13.02564 12.92563 13.12564 12.87270 13.17858
## Mar 2016       13.11409 13.01176 13.21642 12.95758 13.27059
## Apr 2016       13.06203 12.95736 13.16669 12.90196 13.22209
## May 2016       13.03888 12.93197 13.14579 12.87537 13.20239
## Jun 2016       13.00052 12.89139 13.10965 12.83362 13.16743
## Jul 2016       13.07097 12.95968 13.18227 12.90076 13.24119
## Aug 2016       13.07571 12.96229 13.18914 12.90224 13.24919
## Sep 2016       13.02441 12.90889 13.13992 12.84774 13.20107
## Oct 2016       13.06679 12.94922 13.18437 12.88698 13.24661
## Nov 2016       13.04360 12.92401 13.16318 12.86071 13.22648
## Dec 2016       13.17061 13.04903 13.29219 12.98468 13.35655
## Jan 2017       13.05617 12.93330 13.17904 12.86826 13.24408
## Feb 2017       13.03266 12.90822 13.15709 12.84235 13.22296
## Mar 2017       13.12479 12.99868 13.25089 12.93193 13.31764
## Apr 2017       13.08114 12.95346 13.20883 12.88586 13.27642
## May 2017       13.05784 12.92856 13.18712 12.86012 13.25556
## Jun 2017       13.01275 12.88191 13.14360 12.81265 13.21286
## Jul 2017       13.08970 12.95731 13.22210 12.88723 13.29218
## Aug 2017       13.09458 12.96065 13.22850 12.88976 13.29939
## Sep 2017       13.03977 12.90434 13.17521 12.83264 13.24690
## Oct 2017       13.09205 12.95512 13.22898 12.88263 13.30147
## Nov 2017       13.06462 12.92622 13.20303 12.85295 13.27630
## Dec 2017       13.20159 13.06171 13.34147 12.98766 13.41552
## Jan 2018       13.08123 12.93921 13.22324 12.86404 13.29842
## Feb 2018       13.05667 12.91300 13.20033 12.83695 13.27638
## Mar 2018       13.14960 13.00440 13.29479 12.92754 13.37166
## Apr 2018       13.10526 12.95849 13.25203 12.88080 13.32972
## May 2018       13.08215 12.93386 13.23045 12.85536 13.30895
## Jun 2018       13.03535 12.88554 13.18517 12.80623 13.26448
## Jul 2018       13.11345 12.96213 13.26477 12.88202 13.34487
## Aug 2018       13.11880 12.96598 13.27161 12.88509 13.35250
## Sep 2018       13.06455 12.91027 13.21884 12.82860 13.30051
## Oct 2018       13.11613 12.96038 13.27188 12.87794 13.35433
## Nov 2018       13.08906 12.93187 13.24625 12.84866 13.32946
## Dec 2018       13.22767 13.06903 13.38630 12.98505 13.47028
## Jan 2019       13.10596 12.94511 13.26681 12.85996 13.35196
## Feb 2019       13.08224 12.91973 13.24475 12.83371 13.33077
## Mar 2019       13.17465 13.01061 13.33869 12.92378 13.42552
## Apr 2019       13.12917 12.96356 13.29478 12.87589 13.38245
## May 2019       13.10608 12.93894 13.27322 12.85046 13.36170
## Jun 2019       13.06024 12.89157 13.22891 12.80228 13.31820
## Jul 2019       13.13742 12.96724 13.30760 12.87716 13.39768
## Aug 2019       13.14274 12.97106 13.31442 12.88018 13.40530
## Sep 2019       13.08896 12.91581 13.26212 12.82414 13.35379
## Oct 2019       13.13920 12.96456 13.31384 12.87212 13.40628
## Nov 2019       13.11270 12.93661 13.28879 12.84340 13.38200
## Dec 2019       13.24990 13.07235 13.42745 12.97836 13.52144
## Jan 2020       13.12904 12.94947 13.30862 12.85440 13.40368
## Feb 2020       13.10545 12.92428 13.28661 12.82838 13.38251
## Mar 2020       13.19776 13.01510 13.38041 12.91841 13.47710
## Apr 2020       13.15239 12.96821 13.33657 12.87071 13.43407
## May 2020       13.12928 12.94361 13.31495 12.84532 13.41323
## Jun 2020       13.08366 12.89650 13.27081 12.79742 13.36989
## Jul 2020       13.16069 12.97207 13.34932 12.87221 13.44918
## Aug 2020       13.16595 12.97586 13.35604 12.87523 13.45667
## Sep 2020       13.11209 12.92055 13.30363 12.81915 13.40502
## Oct 2020       13.16245 12.96946 13.35543 12.86731 13.45759
## Nov 2020       13.13589 12.94148 13.33029 12.83857 13.43320
## Dec 2020       13.27290 13.07706 13.46873 12.97339 13.57241
## Jan 2021       13.15220 12.95439 13.35001 12.84968 13.45472
## Feb 2021       13.12849 12.92914 13.32784 12.82361 13.43337
## Mar 2021       13.22087 13.02006 13.42168 12.91376 13.52798
## Apr 2021       13.17566 12.97337 13.37796 12.86628 13.48504
## May 2021       13.15255 12.94880 13.35629 12.84094 13.46415

Como se puede ver, el modelo es bueno para predecir tendencias. Sin embargo, no se ve muy acertado para los datos que se tiene en el conjunto de test.

###Prophet algoritmo de FB

## Warning: package 'Rcpp' was built under R version 4.0.5
## Warning: package 'rlang' was built under R version 4.0.5
## Warning: package 'prophet' was built under R version 4.0.5
## Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.

Veamos que tanto se acerca la predicción a los datos reales, en el siguiente gráfico. Los datos del conjunto de prueba se representan en rojo mientras los datos de la predicción se representan en azul.

como podemos observar, las predicciones realizadas por prophet fueron peores y erroneas en comparacion con la funcion de Arima.

##Prediccion de los ultimos 3 años

train <- head(consumoGS_ts, round(length(consumoGS_ts) * 0.891632653))
h <- length(consumoGS_ts) - length(train)
test <- tail(consumoGS_ts, h)
end(train)
## [1] 2019    1
fitArima %>%
  forecast(h) 
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Jan 2015       13.03072 12.96888 13.09256 12.93614 13.12530
## Feb 2015       13.00463 12.93677 13.07249 12.90085 13.10841
## Mar 2015       13.09947 13.02944 13.16950 12.99237 13.20658
## Apr 2015       13.04365 12.96995 13.11735 12.93094 13.15636
## May 2015       13.02189 12.94554 13.09824 12.90512 13.13866
## Jun 2015       12.97006 12.89072 13.04940 12.84872 13.09139
## Jul 2015       13.04983 12.96784 13.13183 12.92443 13.17523
## Aug 2015       13.05805 12.97337 13.14273 12.92854 13.18755
## Sep 2015       13.01030 12.92308 13.09752 12.87690 13.14370
## Oct 2015       13.04912 12.95939 13.13885 12.91189 13.18635
## Nov 2015       13.02783 12.93569 13.11997 12.88691 13.16875
## Dec 2015       13.16829 13.07378 13.26281 13.02374 13.31284
## Jan 2016       13.04320 12.94569 13.14071 12.89407 13.19233
## Feb 2016       13.02564 12.92563 13.12564 12.87270 13.17858
## Mar 2016       13.11409 13.01176 13.21642 12.95758 13.27059
## Apr 2016       13.06203 12.95736 13.16669 12.90196 13.22209
## May 2016       13.03888 12.93197 13.14579 12.87537 13.20239
## Jun 2016       13.00052 12.89139 13.10965 12.83362 13.16743
## Jul 2016       13.07097 12.95968 13.18227 12.90076 13.24119
## Aug 2016       13.07571 12.96229 13.18914 12.90224 13.24919
## Sep 2016       13.02441 12.90889 13.13992 12.84774 13.20107
## Oct 2016       13.06679 12.94922 13.18437 12.88698 13.24661
## Nov 2016       13.04360 12.92401 13.16318 12.86071 13.22648
## Dec 2016       13.17061 13.04903 13.29219 12.98468 13.35655
## Jan 2017       13.05617 12.93330 13.17904 12.86826 13.24408
## Feb 2017       13.03266 12.90822 13.15709 12.84235 13.22296
## Mar 2017       13.12479 12.99868 13.25089 12.93193 13.31764
## Apr 2017       13.08114 12.95346 13.20883 12.88586 13.27642
plot(forecast(fitArima,h))

df<-data.frame(ds=as.Date(as.yearmon(time(train))),y=as.matrix(train) )
testdf<-data.frame(ds=as.Date(as.yearmon(time(test))),y=as.matrix(test) )
fitProphet<-prophet(df,yearly.seasonality = T,weekly.seasonality = T)
## Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.
future <- make_future_dataframe(fitProphet,periods = h,freq = "month", include_history = T)
p <- predict(fitProphet,future)
p<-p[,c("ds","yhat","yhat_lower","yhat_upper")]
plot(fitProphet,p)

pred<-tail(p,h)
pred$y<-testdf$y


ggplot(pred, aes(x=ds, y=yhat)) +
   geom_line(size=1, alpha=0.8) +
   geom_ribbon(aes(ymin=yhat_lower, ymax=yhat_upper), fill="blue", alpha=0.2) +
   geom_line(data=pred, aes(x=ds, y=y),color="red")

En la grafica de los ultimos tres anios podemos ver que fue accurate su predccion con un bajon de consumo de gasolina a mediados de 2020. La prediccion dio lugar a que la tendencia fuera media y nomral, pero en los datos normal (en rojo) podemos ver que se tuvo un decrecimeinto significativo cuando se tuvo la pandemia, lo cual tiene sentido ya que las restricciones no permitian la movilizacion por carro a inicios de esta.

Serie de tiempo Gasolina Super Importacion

Union de columnas de tiempo

Creando la serie de tiempo. Viendo incio, fin y frencuencia

## [1] 2001    1
## [1] 2021    5
## [1] 12

Creando el grafico de la serie

Se observa que es una serie de tiempo no estacionaria ya que su componente tendencia no es constante, si no que es creciente. Ademas, la varianza tampoco es constante.

Descomponiendo la serie

Se puede observar que la componente tendencia es creciente, la componente estacional si tiene un patron y la varianza no es constante. Se tendra que transformar

Transformacion logaritmica para estacionar serie

logSS <- log(trainS)
lambda <- BoxCox.lambda(trainS)
boxcoxGR<-BoxCox(trainS,lambda)

plot(decompose(logSS))

plot(logSS)

Al hacer la transformacion logaritmica, se observa que la serie se ha estabilizado un poco mas en media. Sin embargo en varianza sigue lejos de estar constante

Ahora se haran las pruebas de DickeyFuller

adfTest(trainS)
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: -0.567
##   P VALUE:
##     0.4339 
## 
## Description:
##  Fri Aug 06 21:43:01 2021 by user: jfdel
unitrootTest(trainS)
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     DF: -0.567
##   P VALUE:
##     t: 0.4704 
##     n: 0.5542 
## 
## Description:
##  Fri Aug 06 21:43:01 2021 by user: jfdel

Con las pruebas de Dickey-Fuller, se observa que p es mayor a 0.05. No se puede rechazar Ho y existen raices unitarias. Aunque la serie se vea mas estable en media, aun no es estacionaria en media

Ahora se hace con una diferenciacion en la serie de tiempo

GS_diff<-diff(logSS)
adfTest(diff(trainS))
## Warning in adfTest(diff(trainS)): p-value smaller than printed p-value
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: -16.9539
##   P VALUE:
##     0.01 
## 
## Description:
##  Fri Aug 06 21:43:01 2021 by user: jfdel

Ahora el valor p si es menor a 0.05. Ahora si se puede rechazar la Ho, y ahora no hay raices unitarias

Valores de p y q

Se observa que el ACF se anula muy rapido y cerca de uno por lo que p puede ser uno. En el PACF, se observa que q es 1 ya que hay un pico en el retardo que se sale en 1. Por lo tanto tenemos un modelo ARIMA(1,1,1)

Valores p y q del la componente estacional

p es 1 ya que se tiene solamente un ciclo antes de que se anule q es 1 ya que se acerca a la frecuencia de 12, antes de anularse Por lo tanto, tenemos un modelo ARIMA (1,1,1)

Ahora se comprobara mediante fitarima si en efecto, los coeficiente son significativos

## 
## z test of coefficients:
## 
##       Estimate Std. Error  z value  Pr(>|z|)    
## ar1  -0.418542   0.077040  -5.4328 5.549e-08 ***
## ma1  -0.856811   0.045973 -18.6372 < 2.2e-16 ***
## sar1 -0.197184   0.080253  -2.4570   0.01401 *  
## sma1 -0.999990   0.099506 -10.0495 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## z test of coefficients:
## 
##          Estimate  Std. Error z value  Pr(>|z|)    
## ar1     -0.654059    0.148938 -4.3915 1.126e-05 ***
## ma1     -0.537645    0.181707 -2.9589  0.003088 ** 
## ma2     -0.299850    0.177230 -1.6919  0.090670 .  
## sma1    -0.161559    0.081102 -1.9921  0.046365 *  
## drift 1299.928853  616.505243  2.1085  0.034984 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Como se puede observar, los coeficientes son significativos.

Ahora se hace el analisis de los residuales

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,1)(1,1,1)[12]
## Q* = 19.747, df = 20, p-value = 0.4738
## 
## Model df: 4.   Total lags used: 24

Los residuos se comportan de manera normal y continua. Ademas, en la prueba L-jung-Box se tuvo un valor de p mayor a 0.05. Por lo que es un modelo apto para predecir.

Ahora se hace la prediccion arima en si.

##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## May 2015       13.15250 12.85453 13.45046 12.69680 13.60819
## Jun 2015       13.19007 12.88110 13.49905 12.71754 13.66261
## Jul 2015       13.05966 12.74122 13.37810 12.57264 13.54667
## Aug 2015       13.11468 12.79608 13.43328 12.62742 13.60194
## Sep 2015       13.04087 12.71997 13.36177 12.55010 13.53165
## Oct 2015       13.20716 12.88517 13.52914 12.71472 13.69959
## Nov 2015       13.10977 12.78623 13.43331 12.61496 13.60458
## Dec 2015       13.23458 12.90986 13.55929 12.73797 13.73118
## Jan 2016       13.17263 12.84683 13.49842 12.67437 13.67089
## Feb 2016       13.15992 12.83276 13.48709 12.65957 13.66028
## Mar 2016       13.23094 12.90240 13.55948 12.72848 13.73340
## Apr 2016       13.12845 12.79855 13.45836 12.62391 13.63300
## May 2016       13.24550 12.91553 13.57546 12.74086 13.75013
## Jun 2016       13.17113 12.83867 13.50360 12.66267 13.67960
## Jul 2016       13.14640 12.81334 13.47947 12.63702 13.65578
## Aug 2016       13.11720 12.78290 13.45149 12.60594 13.62846
## Sep 2016       13.12343 12.78819 13.45866 12.61073 13.63613
## Oct 2016       13.18349 12.84722 13.51976 12.66921 13.69777
## Nov 2016       13.17294 12.83565 13.51024 12.65709 13.68879
## Dec 2016       13.29044 12.95228 13.62861 12.77326 13.80762
## Jan 2017       13.19911 12.86005 13.53818 12.68056 13.71766
## Feb 2017       13.25628 12.91622 13.59633 12.73621 13.77635
## Mar 2017       13.30061 12.95957 13.64165 12.77904 13.82219
## Apr 2017       13.21080 12.86878 13.55283 12.68772 13.73388
## May 2017       13.27490 12.92861 13.62119 12.74530 13.80451
## Jun 2017       13.22261 12.87586 13.56936 12.69230 13.75292
## Jul 2017       13.17704 12.82870 13.52538 12.64430 13.70978
## Aug 2017       13.16445 12.81507 13.51382 12.63013 13.69877
## Sep 2017       13.15489 12.80426 13.50552 12.61865 13.69113
## Oct 2017       13.23590 12.88414 13.58766 12.69793 13.77386
## Nov 2017       13.20823 12.85524 13.56121 12.66839 13.74807
## Dec 2017       13.32717 12.97320 13.68114 12.78582 13.86852
## Jan 2018       13.24163 12.88654 13.59673 12.69856 13.78471
## Feb 2018       13.28502 12.92878 13.64126 12.74020 13.82984
## Mar 2018       13.33462 12.97724 13.69200 12.78806 13.88118
## Apr 2018       13.24231 12.88380 13.60082 12.69401 13.79060
## May 2018       13.31685 12.95568 13.67801 12.76449 13.86920
## Jun 2018       13.26020 12.89823 13.62217 12.70662 13.81379
## Jul 2018       13.21874 12.85532 13.58216 12.66294 13.77455
## Aug 2018       13.20287 12.83830 13.56745 12.64531 13.76044
## Sep 2018       13.19643 12.83058 13.56228 12.63691 13.75595
## Oct 2018       13.27331 12.90627 13.64034 12.71197 13.83464
## Nov 2018       13.24901 12.88069 13.61734 12.68571 13.81231
## Dec 2018       13.36767 12.99834 13.73701 12.80282 13.93252
## Jan 2019       13.28099 12.91040 13.65159 12.71421 13.84778
## Feb 2019       13.32710 12.95532 13.69888 12.75851 13.89569
## Mar 2019       13.37566 13.00269 13.74862 12.80526 13.94605
## Apr 2019       13.28384 12.90970 13.65798 12.71164 13.85604
## May 2019       13.35632 12.97916 13.73348 12.77950 13.93314
## Jun 2019       13.30054 12.92258 13.67849 12.72251 13.87857
## Jul 2019       13.25826 12.87875 13.63777 12.67785 13.83868
## Aug 2019       13.24304 12.86233 13.62375 12.66080 13.82528
## Sep 2019       13.23598 12.85392 13.61804 12.65167 13.82030
## Oct 2019       13.31367 12.93037 13.69698 12.72746 13.89989
## Nov 2019       13.28872 12.90405 13.67338 12.70042 13.87701
## Dec 2019       13.40743 13.02171 13.79314 12.81753 13.99733
## Jan 2020       13.32098 12.93386 13.70809 12.72893 13.91302
## Feb 2020       13.36655 12.97819 13.75490 12.77261 13.96048
## Mar 2020       13.41531 13.02572 13.80490 12.81948 14.01114
## Apr 2020       13.32339 12.93257 13.71422 12.72568 13.92111
## May 2020       13.39628 13.00244 13.79013 12.79395 13.99861
## Jun 2020       13.34033 12.94563 13.73502 12.73670 13.94396
## Jul 2020       13.29821 12.90192 13.69451 12.69213 13.90430
## Aug 2020       13.28286 12.88532 13.68041 12.67487 13.89086
## Sep 2020       13.27593 12.87697 13.67488 12.66578 13.88608
## Oct 2020       13.35346 12.95321 13.75370 12.74134 13.96558
## Nov 2020       13.32863 12.92696 13.73030 12.71433 13.94293
## Dec 2020       13.44733 13.04458 13.85008 12.83138 14.06329
## Jan 2021       13.36084 12.95657 13.76511 12.74256 13.97911
## Feb 2021       13.40651 13.00096 13.81207 12.78627 14.02675
## Mar 2021       13.45523 13.04840 13.86207 12.83303 14.07744
## Apr 2021       13.36334 12.95522 13.77145 12.73918 13.98749
## May 2021       13.43615 13.02495 13.84734 12.80728 14.06501

Se puede observar que el algoritmo pudo predecir la tendencia, sin embargo, no pudo predecir los datos.

Analisis con profet de facebook gasolina super importacion

## function (package, help, pos = 2, lib.loc = NULL, character.only = FALSE, 
##     logical.return = FALSE, warn.conflicts, quietly = FALSE, 
##     verbose = getOption("verbose"), mask.ok, exclude, include.only, 
##     attach.required = missing(include.only)) 
## {
##     conf.ctrl <- getOption("conflicts.policy")
##     if (is.character(conf.ctrl)) 
##         conf.ctrl <- switch(conf.ctrl, strict = list(error = TRUE, 
##             warn = FALSE), depends.ok = list(error = TRUE, generics.ok = TRUE, 
##             can.mask = c("base", "methods", "utils", "grDevices", 
##                 "graphics", "stats"), depends.ok = TRUE), warning(gettextf("unknown conflict policy: %s", 
##             sQuote(conf.ctrl)), call. = FALSE, domain = NA))
##     if (!is.list(conf.ctrl)) 
##         conf.ctrl <- NULL
##     stopOnConflict <- isTRUE(conf.ctrl$error)
##     if (missing(warn.conflicts)) 
##         warn.conflicts <- if (isFALSE(conf.ctrl$warn)) 
##             FALSE
##         else TRUE
##     if ((!missing(include.only)) && (!missing(exclude))) 
##         stop(gettext("only one of 'include.only' and 'exclude' can be used"), 
##             call. = FALSE, domain = NA)
##     testRversion <- function(pkgInfo, pkgname, pkgpath) {
##         if (is.null(built <- pkgInfo$Built)) 
##             stop(gettextf("package %s has not been installed properly\n", 
##                 sQuote(pkgname)), call. = FALSE, domain = NA)
##         R_version_built_under <- as.numeric_version(built$R)
##         if (R_version_built_under < "3.0.0") 
##             stop(gettextf("package %s was built before R 3.0.0: please re-install it", 
##                 sQuote(pkgname)), call. = FALSE, domain = NA)
##         current <- getRversion()
##         if (length(Rdeps <- pkgInfo$Rdepends2)) {
##             for (dep in Rdeps) if (length(dep) > 1L) {
##                 target <- dep$version
##                 res <- do.call(dep$op, if (is.character(target)) 
##                   list(as.numeric(R.version[["svn rev"]]), as.numeric(sub("^r", 
##                     "", target)))
##                 else list(current, as.numeric_version(target)))
##                 if (!res) 
##                   stop(gettextf("This is R %s, package %s needs %s %s", 
##                     current, sQuote(pkgname), dep$op, target), 
##                     call. = FALSE, domain = NA)
##             }
##         }
##         if (R_version_built_under > current) 
##             warning(gettextf("package %s was built under R version %s", 
##                 sQuote(pkgname), as.character(built$R)), call. = FALSE, 
##                 domain = NA)
##         platform <- built$Platform
##         r_arch <- .Platform$r_arch
##         if (.Platform$OS.type == "unix") {
##         }
##         else {
##             if (nzchar(platform) && !grepl("mingw", platform)) 
##                 stop(gettextf("package %s was built for %s", 
##                   sQuote(pkgname), platform), call. = FALSE, 
##                   domain = NA)
##         }
##         if (nzchar(r_arch) && file.exists(file.path(pkgpath, 
##             "libs")) && !file.exists(file.path(pkgpath, "libs", 
##             r_arch))) 
##             stop(gettextf("package %s is not installed for 'arch = %s'", 
##                 sQuote(pkgname), r_arch), call. = FALSE, domain = NA)
##     }
##     checkNoGenerics <- function(env, pkg) {
##         nenv <- env
##         ns <- .getNamespace(as.name(pkg))
##         if (!is.null(ns)) 
##             nenv <- asNamespace(ns)
##         if (exists(".noGenerics", envir = nenv, inherits = FALSE)) 
##             TRUE
##         else {
##             !any(startsWith(names(env), ".__T"))
##         }
##     }
##     checkConflicts <- function(package, pkgname, pkgpath, nogenerics, 
##         env) {
##         dont.mind <- c("last.dump", "last.warning", ".Last.value", 
##             ".Random.seed", ".Last.lib", ".onDetach", ".packageName", 
##             ".noGenerics", ".required", ".no_S3_generics", ".Depends", 
##             ".requireCachedGenerics")
##         sp <- search()
##         lib.pos <- which(sp == pkgname)
##         ob <- names(as.environment(lib.pos))
##         if (!nogenerics) {
##             these <- ob[startsWith(ob, ".__T__")]
##             gen <- gsub(".__T__(.*):([^:]+)", "\\1", these)
##             from <- gsub(".__T__(.*):([^:]+)", "\\2", these)
##             gen <- gen[from != package]
##             ob <- ob[!(ob %in% gen)]
##         }
##         ipos <- seq_along(sp)[-c(lib.pos, match(c("Autoloads", 
##             "CheckExEnv"), sp, 0L))]
##         cpos <- NULL
##         conflicts <- vector("list", 0)
##         for (i in ipos) {
##             obj.same <- match(names(as.environment(i)), ob, nomatch = 0L)
##             if (any(obj.same > 0L)) {
##                 same <- ob[obj.same]
##                 same <- same[!(same %in% dont.mind)]
##                 Classobjs <- which(startsWith(same, ".__"))
##                 if (length(Classobjs)) 
##                   same <- same[-Classobjs]
##                 same.isFn <- function(where) vapply(same, exists, 
##                   NA, where = where, mode = "function", inherits = FALSE)
##                 same <- same[same.isFn(i) == same.isFn(lib.pos)]
##                 not.Ident <- function(ch, TRAFO = identity, ...) vapply(ch, 
##                   function(.) !identical(TRAFO(get(., i)), TRAFO(get(., 
##                     lib.pos)), ...), NA)
##                 if (length(same)) 
##                   same <- same[not.Ident(same)]
##                 if (length(same) && identical(sp[i], "package:base")) 
##                   same <- same[not.Ident(same, ignore.environment = TRUE)]
##                 if (length(same)) {
##                   conflicts[[sp[i]]] <- same
##                   cpos[sp[i]] <- i
##                 }
##             }
##         }
##         if (length(conflicts)) {
##             if (stopOnConflict) {
##                 emsg <- ""
##                 pkg <- names(conflicts)
##                 notOK <- vector("list", 0)
##                 for (i in seq_along(conflicts)) {
##                   pkgname <- sub("^package:", "", pkg[i])
##                   if (pkgname %in% canMaskEnv$canMask) 
##                     next
##                   same <- conflicts[[i]]
##                   if (is.list(mask.ok)) 
##                     myMaskOK <- mask.ok[[pkgname]]
##                   else myMaskOK <- mask.ok
##                   if (isTRUE(myMaskOK)) 
##                     same <- NULL
##                   else if (is.character(myMaskOK)) 
##                     same <- setdiff(same, myMaskOK)
##                   if (length(same)) {
##                     notOK[[pkg[i]]] <- same
##                     msg <- .maskedMsg(sort(same), pkg = sQuote(pkg[i]), 
##                       by = cpos[i] < lib.pos)
##                     emsg <- paste(emsg, msg, sep = "\n")
##                   }
##                 }
##                 if (length(notOK)) {
##                   msg <- gettextf("Conflicts attaching package %s:\n%s", 
##                     sQuote(package), emsg)
##                   stop(errorCondition(msg, package = package, 
##                     conflicts = conflicts, class = "packageConflictError"))
##                 }
##             }
##             if (warn.conflicts) {
##                 packageStartupMessage(gettextf("\nAttaching package: %s\n", 
##                   sQuote(package)), domain = NA)
##                 pkg <- names(conflicts)
##                 for (i in seq_along(conflicts)) {
##                   msg <- .maskedMsg(sort(conflicts[[i]]), pkg = sQuote(pkg[i]), 
##                     by = cpos[i] < lib.pos)
##                   packageStartupMessage(msg, domain = NA)
##                 }
##             }
##         }
##     }
##     if (verbose && quietly) 
##         message("'verbose' and 'quietly' are both true; being verbose then ..")
##     if (!missing(package)) {
##         if (is.null(lib.loc)) 
##             lib.loc <- .libPaths()
##         lib.loc <- lib.loc[dir.exists(lib.loc)]
##         if (!character.only) 
##             package <- as.character(substitute(package))
##         if (length(package) != 1L) 
##             stop("'package' must be of length 1")
##         if (is.na(package) || (package == "")) 
##             stop("invalid package name")
##         pkgname <- paste0("package:", package)
##         newpackage <- is.na(match(pkgname, search()))
##         if (newpackage) {
##             pkgpath <- find.package(package, lib.loc, quiet = TRUE, 
##                 verbose = verbose)
##             if (length(pkgpath) == 0L) {
##                 if (length(lib.loc) && !logical.return) 
##                   stop(packageNotFoundError(package, lib.loc, 
##                     sys.call()))
##                 txt <- if (length(lib.loc)) 
##                   gettextf("there is no package called %s", sQuote(package))
##                 else gettext("no library trees found in 'lib.loc'")
##                 if (logical.return) {
##                   warning(txt, domain = NA)
##                   return(FALSE)
##                 }
##                 else stop(txt, domain = NA)
##             }
##             which.lib.loc <- normalizePath(dirname(pkgpath), 
##                 "/", TRUE)
##             pfile <- system.file("Meta", "package.rds", package = package, 
##                 lib.loc = which.lib.loc)
##             if (!nzchar(pfile)) 
##                 stop(gettextf("%s is not a valid installed package", 
##                   sQuote(package)), domain = NA)
##             pkgInfo <- readRDS(pfile)
##             testRversion(pkgInfo, package, pkgpath)
##             if (is.character(pos)) {
##                 npos <- match(pos, search())
##                 if (is.na(npos)) {
##                   warning(gettextf("%s not found on search path, using pos = 2", 
##                     sQuote(pos)), domain = NA)
##                   pos <- 2
##                 }
##                 else pos <- npos
##             }
##             deps <- unique(names(pkgInfo$Depends))
##             depsOK <- isTRUE(conf.ctrl$depends.ok)
##             if (depsOK) {
##                 canMaskEnv <- dynGet("__library_can_mask__", 
##                   NULL)
##                 if (is.null(canMaskEnv)) {
##                   canMaskEnv <- new.env()
##                   canMaskEnv$canMask <- union("base", conf.ctrl$can.mask)
##                   "__library_can_mask__" <- canMaskEnv
##                 }
##                 canMaskEnv$canMask <- unique(c(package, deps, 
##                   canMaskEnv$canMask))
##             }
##             else canMaskEnv <- NULL
##             if (attach.required) 
##                 .getRequiredPackages2(pkgInfo, quietly = quietly)
##             cr <- conflictRules(package)
##             if (missing(mask.ok)) 
##                 mask.ok <- cr$mask.ok
##             if (missing(exclude)) 
##                 exclude <- cr$exclude
##             if (packageHasNamespace(package, which.lib.loc)) {
##                 if (isNamespaceLoaded(package)) {
##                   newversion <- as.numeric_version(pkgInfo$DESCRIPTION["Version"])
##                   oldversion <- as.numeric_version(getNamespaceVersion(package))
##                   if (newversion != oldversion) {
##                     tryCatch(unloadNamespace(package), error = function(e) {
##                       P <- if (!is.null(cc <- conditionCall(e))) 
##                         paste("Error in", deparse(cc)[1L], ": ")
##                       else "Error : "
##                       stop(gettextf("Package %s version %s cannot be unloaded:\n %s", 
##                         sQuote(package), oldversion, paste0(P, 
##                           conditionMessage(e), "\n")), domain = NA)
##                     })
##                   }
##                 }
##                 tt <- tryCatch({
##                   attr(package, "LibPath") <- which.lib.loc
##                   ns <- loadNamespace(package, lib.loc)
##                   env <- attachNamespace(ns, pos = pos, deps, 
##                     exclude, include.only)
##                 }, error = function(e) {
##                   P <- if (!is.null(cc <- conditionCall(e))) 
##                     paste(" in", deparse(cc)[1L])
##                   else ""
##                   msg <- gettextf("package or namespace load failed for %s%s:\n %s", 
##                     sQuote(package), P, conditionMessage(e))
##                   if (logical.return) 
##                     message(paste("Error:", msg), domain = NA)
##                   else stop(msg, call. = FALSE, domain = NA)
##                 })
##                 if (logical.return && is.null(tt)) 
##                   return(FALSE)
##                 attr(package, "LibPath") <- NULL
##                 {
##                   on.exit(detach(pos = pos))
##                   nogenerics <- !.isMethodsDispatchOn() || checkNoGenerics(env, 
##                     package)
##                   if (isFALSE(conf.ctrl$generics.ok) || (stopOnConflict && 
##                     !isTRUE(conf.ctrl$generics.ok))) 
##                     nogenerics <- TRUE
##                   if (stopOnConflict || (warn.conflicts && !exists(".conflicts.OK", 
##                     envir = env, inherits = FALSE))) 
##                     checkConflicts(package, pkgname, pkgpath, 
##                       nogenerics, ns)
##                   on.exit()
##                   if (logical.return) 
##                     return(TRUE)
##                   else return(invisible(.packages()))
##                 }
##             }
##             else stop(gettextf("package %s does not have a namespace and should be re-installed", 
##                 sQuote(package)), domain = NA)
##         }
##         if (verbose && !newpackage) 
##             warning(gettextf("package %s already present in search()", 
##                 sQuote(package)), domain = NA)
##     }
##     else if (!missing(help)) {
##         if (!character.only) 
##             help <- as.character(substitute(help))
##         pkgName <- help[1L]
##         pkgPath <- find.package(pkgName, lib.loc, verbose = verbose)
##         docFiles <- c(file.path(pkgPath, "Meta", "package.rds"), 
##             file.path(pkgPath, "INDEX"))
##         if (file.exists(vignetteIndexRDS <- file.path(pkgPath, 
##             "Meta", "vignette.rds"))) 
##             docFiles <- c(docFiles, vignetteIndexRDS)
##         pkgInfo <- vector("list", 3L)
##         readDocFile <- function(f) {
##             if (basename(f) %in% "package.rds") {
##                 txt <- readRDS(f)$DESCRIPTION
##                 if ("Encoding" %in% names(txt)) {
##                   to <- if (Sys.getlocale("LC_CTYPE") == "C") 
##                     "ASCII//TRANSLIT"
##                   else ""
##                   tmp <- try(iconv(txt, from = txt["Encoding"], 
##                     to = to))
##                   if (!inherits(tmp, "try-error")) 
##                     txt <- tmp
##                   else warning("'DESCRIPTION' has an 'Encoding' field and re-encoding is not possible", 
##                     call. = FALSE)
##                 }
##                 nm <- paste0(names(txt), ":")
##                 formatDL(nm, txt, indent = max(nchar(nm, "w")) + 
##                   3L)
##             }
##             else if (basename(f) %in% "vignette.rds") {
##                 txt <- readRDS(f)
##                 if (is.data.frame(txt) && nrow(txt)) 
##                   cbind(basename(gsub("\\.[[:alpha:]]+$", "", 
##                     txt$File)), paste(txt$Title, paste0(rep.int("(source", 
##                     NROW(txt)), ifelse(nzchar(txt$PDF), ", pdf", 
##                     ""), ")")))
##                 else NULL
##             }
##             else readLines(f)
##         }
##         for (i in which(file.exists(docFiles))) pkgInfo[[i]] <- readDocFile(docFiles[i])
##         y <- list(name = pkgName, path = pkgPath, info = pkgInfo)
##         class(y) <- "packageInfo"
##         return(y)
##     }
##     else {
##         if (is.null(lib.loc)) 
##             lib.loc <- .libPaths()
##         db <- matrix(character(), nrow = 0L, ncol = 3L)
##         nopkgs <- character()
##         for (lib in lib.loc) {
##             a <- .packages(all.available = TRUE, lib.loc = lib)
##             for (i in sort(a)) {
##                 file <- system.file("Meta", "package.rds", package = i, 
##                   lib.loc = lib)
##                 title <- if (nzchar(file)) {
##                   txt <- readRDS(file)
##                   if (is.list(txt)) 
##                     txt <- txt$DESCRIPTION
##                   if ("Encoding" %in% names(txt)) {
##                     to <- if (Sys.getlocale("LC_CTYPE") == "C") 
##                       "ASCII//TRANSLIT"
##                     else ""
##                     tmp <- try(iconv(txt, txt["Encoding"], to, 
##                       "?"))
##                     if (!inherits(tmp, "try-error")) 
##                       txt <- tmp
##                     else warning("'DESCRIPTION' has an 'Encoding' field and re-encoding is not possible", 
##                       call. = FALSE)
##                   }
##                   txt["Title"]
##                 }
##                 else NA
##                 if (is.na(title)) 
##                   title <- " ** No title available ** "
##                 db <- rbind(db, cbind(i, lib, title))
##             }
##             if (length(a) == 0L) 
##                 nopkgs <- c(nopkgs, lib)
##         }
##         dimnames(db) <- list(NULL, c("Package", "LibPath", "Title"))
##         if (length(nopkgs) && !missing(lib.loc)) {
##             pkglist <- paste(sQuote(nopkgs), collapse = ", ")
##             msg <- sprintf(ngettext(length(nopkgs), "library %s contains no packages", 
##                 "libraries %s contain no packages"), pkglist)
##             warning(msg, domain = NA)
##         }
##         y <- list(header = NULL, results = db, footer = NULL)
##         class(y) <- "libraryIQR"
##         return(y)
##     }
##     if (logical.return) 
##         TRUE
##     else invisible(.packages())
## }
## <bytecode: 0x000000001309d8c8>
## <environment: namespace:base>
## Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.

Como se puede observar, el ARIMA tuvo una mejor prediccion.

Prediccion ultimos 3 anios

## [1] 2019    2
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## May 2015       13.15250 12.85453 13.45046 12.69680 13.60819
## Jun 2015       13.19007 12.88110 13.49905 12.71754 13.66261
## Jul 2015       13.05966 12.74122 13.37810 12.57264 13.54667
## Aug 2015       13.11468 12.79608 13.43328 12.62742 13.60194
## Sep 2015       13.04087 12.71997 13.36177 12.55010 13.53165
## Oct 2015       13.20716 12.88517 13.52914 12.71472 13.69959
## Nov 2015       13.10977 12.78623 13.43331 12.61496 13.60458
## Dec 2015       13.23458 12.90986 13.55929 12.73797 13.73118
## Jan 2016       13.17263 12.84683 13.49842 12.67437 13.67089
## Feb 2016       13.15992 12.83276 13.48709 12.65957 13.66028
## Mar 2016       13.23094 12.90240 13.55948 12.72848 13.73340
## Apr 2016       13.12845 12.79855 13.45836 12.62391 13.63300
## May 2016       13.24550 12.91553 13.57546 12.74086 13.75013
## Jun 2016       13.17113 12.83867 13.50360 12.66267 13.67960
## Jul 2016       13.14640 12.81334 13.47947 12.63702 13.65578
## Aug 2016       13.11720 12.78290 13.45149 12.60594 13.62846
## Sep 2016       13.12343 12.78819 13.45866 12.61073 13.63613
## Oct 2016       13.18349 12.84722 13.51976 12.66921 13.69777
## Nov 2016       13.17294 12.83565 13.51024 12.65709 13.68879
## Dec 2016       13.29044 12.95228 13.62861 12.77326 13.80762
## Jan 2017       13.19911 12.86005 13.53818 12.68056 13.71766
## Feb 2017       13.25628 12.91622 13.59633 12.73621 13.77635
## Mar 2017       13.30061 12.95957 13.64165 12.77904 13.82219
## Apr 2017       13.21080 12.86878 13.55283 12.68772 13.73388
## May 2017       13.27490 12.92861 13.62119 12.74530 13.80451
## Jun 2017       13.22261 12.87586 13.56936 12.69230 13.75292
## Jul 2017       13.17704 12.82870 13.52538 12.64430 13.70978

## Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.

Se puede observar que la prediccion para el 2021 no es aceptable. Ademas, por el efecto de la pandemia y su comportamiento irregular, la prediccion tuvo peores resultados despues de la misma.

####Linea de tiempo para consumo de gasolina Regular Importacion

consumoGR_ts <- xts(consumoRegular$cantRegular, consumoRegular$regularDate) 
consumoGR_ts <-ts_ts(consumoGR_ts)
#Saber cuando empieza la serie y cuando termina
start(consumoGR_ts)
## [1] 2000    1
end(consumoGR_ts)
## [1] 2021    5
#Saber la frecuencia de la serie
frequency(consumoGR_ts)
## [1] 12
#Ver el gráfico de la serie
plot(consumoGR_ts)

###Analisis Importacion Gasolina Regular y Diesel
datosImport <- read.csv(file = 'DatosImportacionCombustibles.csv')

datosImport$Dia <- 1
datosImport$Date <- as.Date(with(datosImport, paste(Anio, Mes, Dia,sep="-")), "%Y-%m-%d")

cantRegular <- c(datosImport$GasolinaRegular)
regularDate <- c(datosImport$Date)
importRegular <- data.frame(cantRegular,regularDate)

cantDiesel <- c(datosImport$TotalDiesel)
dieselDate <- c(datosImport$Date)
importDiesel <- data.frame(cantDiesel,dieselDate)

Inicio y fin de la importacion Gasolina Regular

importRegular_ts <- xts(importRegular$cantRegular, importRegular$regularDate) 
importRegular_ts <-ts_ts(importRegular_ts)

start(importRegular_ts)
## [1] 2001    1
end(importRegular_ts) 
## [1] 2021    5
train <- head(importRegular_ts, round(length(importRegular_ts) * 0.7))
h <- length(importRegular_ts) - length(train)
test <- tail(importRegular_ts, h)

Frecuencia serie importacion Gasolina Regular

frequency(importRegular_ts)
## [1] 12

Grafico de la serie importacion Gasolina Regular

plot(importRegular_ts)
abline(reg=lm(importRegular_ts~time(importRegular_ts)), col=c("red"))

De la grafica puedo observar que tiene tendencia a aumentar. No es estacionaria media

Descomposicion serie importacion Gasolina Regular

dec.GS<-decompose(importRegular_ts)
plot(dec.GS)

Por la descomposicion tampoco parece ser estacionaria en varianza. El factor aleatorio parece tener un patron por lo que no sera una serie util para predecir. Es necesario transformarla

Transformacion log y grafico autocorrelacion serie importacion Gasolina Regular

logIR <- log(train)
lambda <- BoxCox.lambda(train)
boxcoxGR<-BoxCox(train,lambda)

plot(decompose(logIR))

plot(logIR)

acf(logIR)

No hay estacionariedad media

Prueba Dickey-Fuller

adfTest(train)
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: -0.6179
##   P VALUE:
##     0.4177 
## 
## Description:
##  Fri Aug 06 21:43:06 2021 by user: jfdel
unitrootTest(train)
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     DF: -0.6179
##   P VALUE:
##     t: 0.4484 
##     n: 0.5437 
## 
## Description:
##  Fri Aug 06 21:43:06 2021 by user: jfdel

El test Dickey-Fuller resulta con valores mayores a 0.05 por lo que Ho no se rechaza. Si hay raices unitarias. Se debe hacer una diferenciacion de la serie de tiempo.

Diferenciacion serie importacion Gasolina Regular

IR_diff<-diff(logIR)
adfTest(diff(train))
## Warning in adfTest(diff(train)): p-value smaller than printed p-value
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: -15.4037
##   P VALUE:
##     0.01 
## 
## Description:
##  Fri Aug 06 21:43:06 2021 by user: jfdel

Al hacer este cambio el valor p si es menor a 0.05, se rechaza Ho y ahora no hay raices unitarias

Valores p y q

acf(logIR,50) 

pacf(logIR,50) 

Inicio y fin de la importacion Gasolina Diesel

importDiesel_ts <- xts(importDiesel$cantDiesel, importDiesel$dieselDate) 
importDiesel_ts <-ts_ts(importDiesel_ts)

start(importDiesel_ts)
## [1] 2001    1
end(importDiesel_ts) 
## [1] 2021    5
train <- head(importDiesel_ts, round(length(importDiesel_ts) * 0.7))
h <- length(importDiesel_ts) - length(train)
test <- tail(importDiesel_ts, h)

Frecuencia serie importacion Gasolina Diesel

frequency(importDiesel_ts)
## [1] 12

Grafico de la serie importacion Gasolina Diesel

plot(importDiesel_ts)
abline(reg=lm(importDiesel_ts~time(importDiesel_ts)), col=c("red"))

De la grafica puedo observar que tiene tendencia a aumentar. No es estacionaria media

Descomposicion serie importacion Gasolina Diesel

dec.GD<-decompose(importDiesel_ts)
plot(dec.GD)

Por la descomposicion tampoco parece ser estacionaria en varianza. El factor aleatorio parece tener un patron por lo que no sera una serie util para predecir. Es necesario transformarla

Transformacion log y grafico autocorrelacion serie importacion Gasolina Diesel

logID <- log(train)
lambda <- BoxCox.lambda(train)
boxcoxGD<-BoxCox(train,lambda)

plot(decompose(logID))

plot(logID)

acf(logID)

No hay estacionariedad media

Prueba Dickey-Fuller

adfTest(train)
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: -0.9005
##   P VALUE:
##     0.3276 
## 
## Description:
##  Fri Aug 06 21:43:06 2021 by user: jfdel
unitrootTest(train)
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     DF: -0.9005
##   P VALUE:
##     t: 0.3248 
##     n: 0.4892 
## 
## Description:
##  Fri Aug 06 21:43:06 2021 by user: jfdel

El test Dickey-Fuller resulta con valores mayores a 0.05 por lo que Ho se rechaza. No hay raices unitarias. Se debe hacer una diferenciacion de la serie de tiempo.

Diferenciacion serie importacion Gasolina Diesel

ID_diff<-diff(logID)
adfTest(diff(train))
## Warning in adfTest(diff(train)): p-value smaller than printed p-value
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: -15.2594
##   P VALUE:
##     0.01 
## 
## Description:
##  Fri Aug 06 21:43:06 2021 by user: jfdel

Al hacer este cambio el valor p si es menor a 0.05, ya se tienen raices unitarias.

Valores p y q

acf(logID,50) 

pacf(logID,50) 

Analisis Consumo

## Warning: package 'rstan' was built under R version 4.0.5
## Loading required package: StanHeaders
## Warning: package 'StanHeaders' was built under R version 4.0.5
## rstan (Version 2.21.2, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## Do not specify '-march=native' in 'LOCAL_CPPFLAGS' or a Makevars file
## 
## Attaching package: 'rstan'
## The following object is masked from 'package:epiDisplay':
## 
##     lookup

Análisis exploratorio - Box-Plot Diesel, Regular

Luego de realizar el analisis exploratorio de las variables continuas, en relacion a los tipos de gasolina Diesel, Regular y Super, se observo que estas tienen una distribucion normal en base a los diagramas de caja mostrados anteriormente. En ciertos diagramas se puede apreciar sesgos positivos o negativos, sin embargo, estos no son representativos para indicar que no muestran una distribucion normal. En base a lo anterior, se puede indicar que se puede obtener una serie de tiempo del conjunto de datos. Es necesario determinar si la serie de tiempo del conjunto de datos presenta estacionalidad, asi como tambien determinar si se debe realizar una normalizacion en los datos para poder brindar una mejor prediccion que se lograra con los modelos.

Series de Tiempo

Al encontrar las series de tiempo se obtuvieron los datos de las graficas anteriores. Sin embargo es necesario descomponer estas:

Al descomponer la serie de tiempo se puede observar que ninguna de las series son estacionarias ya que la media no es constante al igual que la varianza.

Estacionariedad

Para deterimnar la estacionariedad se realizaron los siguientes graficos que muestra la tendencia de cada serie, factor aleatorio y estacionariedad.

En las graficas anteriores se puede apreicar como estas siguen un patron, esto es indicador de que cumplen con estacionariedad. Se podria indicar que la serie de tiempo del diesel es estacionaria en media, ya que se mantiene relativamente constante, sin embargo, en los ultimos años esta oresenta una baja. Mientras que la serie de tiempo de la gasolina regular no es estacionaria en media.

## Warning in adf.test(logDieselST): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  logDieselST
## Dickey-Fuller = -5.6671, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary
## 
##  Augmented Dickey-Fuller Test
## 
## data:  logRegularST
## Dickey-Fuller = -1.8568, Lag order = 6, p-value = 0.6361
## alternative hypothesis: stationary

Modelos

Modelo Diesel

Modelo Regular

Prophet

Gasolina diesel

## Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.

Gasolina Regular

##           ds        y
## 1 2001-01-01 202645.2
## 2 2001-02-01 205531.0
## 3 2001-03-01 229499.6
## 4 2001-04-01 210680.4
## 5 2001-05-01 208164.3
## 6 2001-06-01 195088.7
## Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.

En base a las graficas anteriores se puede indicar que el modelo fue eficiente Dado que se ajusto el conjunto de datos lo mas posible, se puede observar que hay una cercania entre la prediccion y los datos reales. El mejor modelo producido es el de diesel ya que tiene menos variacion, mientras que la gasolina regular se dificulto la prediccion.